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Abstract 

We have developed a Scalable Linear Augmented Slater- Type Orbital 
(LASTO) method for electronic-structure calculations on free-standing 
atomic clusters. As with other linear methods we solve the Schrodinger 
equation using a mixed basis set consisting of numerical functions inside 
atom-centered spheres and matched onto tail functions outside. The tail 
functions are Slater-type orbitals, which are localized, exponentially de- 
caying functions. To solve the Poisson equation between spheres, we use a 
finite difference method replacing the rapidly varying charge density inside 
the spheres with a smoothed density with the same multipole moments. 
We use multigrid techniques on the mesh, which yield the Coulomb poten- 
tial on the spheres and in turn defines the potential inside via a Dirichlet 
problem. To solve the linear eigen-problem, we use ScaLAPACK, a well- 
developed package to solve large eigensystems with dense matrices. We 
have tested the method on small clusters of palladium. 

1 Introduction 

There are many ways to solve the coupled Schrodinger and Poisson equations 
required for density functional theory. They generally fall into two classes: (a) 
finite cluster calculations using basis sets, such as Gaussians, or (b) periodic 
crystal calculations, as for example those which use a plane wave basis set. 
Within either class the nuclear attraction may be replaced by a pseudo or effec- 
tive core potential. 
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The ability to fabricate nanoscale clusters with tens of thousands of atoms 
has driven a renewed interest in electronic structure methods capable of reaching 
this size. In our research we have opted for the finite cluster approach, since 
many of the systems of interest fall into this category However, the standard 
techniques using Gaussians often require extremely large basis sets and are 
difficult to apply to heavy atoms. Methods for treating crystals which derive 
from the augmented plane wave (APW) method are known to be extremely 
accurate and suitable for all atoms in the periodic table but make essential use 
of Fourier decompositions in a way that does not scale well to massively parallel 
machines. In addition, they make use of supercells, which may not describe 
finite systems accurately. 

Some years ago we developed an augmented basis set method which uses 
linearized solutions of the Schrodinger (or Dirac) equation inside atom centered 
spheres, and Slater-type-orbitals in the region between spheres [UEIE]. It was 
based on Anderson's linear muffin tin orbital method (LMTO) [4] but used 
Slater type orbitals as tail functions in place of Anderson's Spherical Bessel 
functions. While the method was formulated in real space, it was more cleanly 
coded in reciprocal space, which for crystals with small unit cells was equally 
efficient. 

We return here to the real space formulation and apply it to free standing 
clusters. There are a number of other real space formulations of density func- 
tional theory (for a review see [5]). Most of these solve the Schrodinger equation 
directly on a mesh, as for example Chelikowski et al [7] . 

An essential difference is the method for solving the Poisson equation. The 
standard method in crystals [3 [9] is to replace the charge density inside the 
spheres by a smooth pseudo-density which has the same multipole moments as 
the true density (including the nuclear charge) and is represented by a Fourier 
series. Here wc use the same idea, except we represent the pseudo-density on a 
(nonuniform) grid to enhance resolution adaptively, not globally. In addition, we 
use the pseudo-density only for the spherical (monopole) portion of the density, 
as the nonspherical terms can be included directly on the grid. 

In either case, one finds the solution of the Poisson equation outside the 
spheres and interpolates onto the sphere boundaries to define a Dirichlet prob- 
lem for the potential inside which can be solved using the true density. 



2 Methods 

In this section, we consider the Kohn-Sham approach in real space for electronic 
structure calculations. 

We consider the Schrodinger equation 
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V(p(x)) 



Vfc(x) = e k ip k (x), 



(2.1) 
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where 

P (X)= J2 l^(x)| 2 +Pcorc(^(p(x))) (2.2) 



e k <E Fm . 



mi 



is the electronic charge density and 

V(p(x)) = ^Ext(x) + Vkartree(p(x)) + Kcc(p(x)). (2.3) 

is the potential, which consists of the external potential V^Ext, the Hartree po- 
tential Vnartree, and the exchange and correlation potential V xc . We use atomic 
units where lengths are measured in units of the Bohr radius a^H 2 /me 2 = 0.529 
A and energies in Hartrees e 2 /ao = 27.212 eV. The external potential Vjjxt 
is typically a sum of nuclear potentials centered at the atomic positions. The 
Hartree potential Vnartree can be obtained by solving the Poisson equation 



with 



where 



V 2 y H artree(x) = 4^p(x) (2.4) 



lim VHartree(x) = - lim VExt(x), (2.5) 
|x| — >oo |x| — >oc 



V 2 1/ E xt(x) = 4ir Z j S je=0j (x), 
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where Cj's are the center of atoms and Zj's are the nuclear charge. We define 
the Coulomb potential Vc(x) = VExt(x) + Vkartrcc(/0(x)) and solve the Poisson 
equation 

V 2 y c (x) = 4tt(p(x) - ^«x= Cj (x)) (2.6) 

3 

with 

lim Vc(x) = 0. (2.7) 

x| — >oo 

The exchange and correlation potential V xc is formally defined through the 
functional derivative of the exchange and correlation energy, 



SR 

Vyrc = 



Sp 



(2.8) 



In principle, only the density functionals for the exchange and correlation term 
remain to be approximated. Many approximations for the exchange-correlation 
functionals have been developed. We use local spin density (LSD) of the Hedin- 
Lundqvist form [10] and generalized gradient approximations (GGA) of Perdew, 
Burke, Wang and Ernzerhof, (PBE-GGA) [TT] . The GGA approximation is gen- 
erally considered to be more accurate and we confirm this view in our work. 

Due to the functional dependence of V on the density, these equations form a 
set of nonlinear coupled equations. The standard solution procedure is to iterate 
on the solution of the linear subsystems until self-consistency is achieved. We 
show the basic flow-chart to solve the Schrodinger equation with DFT in Fig. 
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Figure 1: Flow chart for the solving Schrodinger equation with DFT * 



To solve the Schrodinger equation, we use the linear augmcntcd-Slater-type- 
orbital (LASTO) basis set which has been described in detail in [IJ[2l[3]. Here, 
we summarize results for the LASTO basis set for numerical implementation. 

We introduce a sphere Sj around the jth atom (Fig. [2]). We use different 
functional forms for the basis functions, 4>inim which is defined with reference 
atom i, inside the Sj and the region outside all Sj, called the interstitial region; 



" *e ^' i yimifi)i in the interstitial region 



T,j\,,[Pinlm,j\ fJ ,9jx(rj) + a in i m ,jxngj\(rj)}yjx^(fj), 



inside Sj 



where Tj and fj 



(2.9) 

(0,(j)) are the spherical coordinates for x = (x, y, z) = 
(r sin 6* cos <f>, r sin 9 sin^, r cos 9) with respect to the position Xj of atom j, Sj 
is a sphere centered at the atom with radius Rj , and the yi m are real spherical 
harmonics. 

In (|2.9|) . the gj\ are numerical solutions of the scalar relativistic radial Dirac 
equation and the gj\ are their energy derivatives. They satisfy the radial equa- 
tions 

h r 9jX = Cj9j\ ( 2 -10) 

and 

h r9j\ = e j9 3 x + 9jX, (2-11) 
where h r is the scalar relativistic radial Hamiltonian, 



2AI 



dr 2 



2 d 
r dr 



A(A + 1) 



1 dVn d 



2M2Mc 2 dr dr 



+ V , 



M = m 



2mc 2 
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Figure 2: Simplified domain in 2D. 



and Vo is the spherical average of the potential. The g's are normalized within 
the spheres, 

Rs 

r 2 g 2 j\{r)dr = 1, 



and gj\ and gj\ are orthogonal. 

The /3's and a's are chosen by matching the interior and the exterior func- 
tions and their derivatives on the boundaries of the spheres, i.e., 

0mim(x) = yjPinlm,j\ti9j\{Tj) + a inhn j x^j A (rj )]j/A/i (fj ) 
A / 1 

= r^e-^'yimih) = C Unm,j x»(rj)yx»{rj), 

X/i 

-j— ^inim(x) = Y^inlmJXfig'jxirj) + OtinlmjXvg'jX )}VXfi (fj ) 

i Xfj, 

for rj = Computation of the Cii nm .jx^ and C inlm require the expansion 
of an STO about a site other than the one on which it is centered. This problem 
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has been considered by many authors pQ. Here we only summarize the results: 

J^Er IrQ™, I'm', Z"m")VJ,,„(r, Rji)yv m »{Rji), for j ^ i, 
OwmjxM - | rr i e -fr tttm(fl)i fori = i) 

(2.12) 

/4irEr^(im > /W,rW0V3 , ,?,(r,ii,-i)W" W .''(^i<), f ° r J # *, 
mZmjAM \(n- l)(rr a e-^ - rr'C^Jttm^), for J = i, 

(2.13) 

where In(lm,l'm' ,l"m") is a Gaunt integral 

I R (lm,l'm',l"m") = / yi m yi<m>yi" m "dr 2 , 



(Rji,Rji) are the spherical coordinates for X, — X,-, and VJ™;/, and are 
given by 

yn (rR) = ( n n-l 1 y W'+i+l) ^ 1^" + j + 1) 1 

' J 1 J c™- 1 ^r(i + i)r(i'-i + i) ^r(j + i)r(/"-j + i)2 i +j+ 1 

x [exp(Cr) + (-l) r+fe - fc '- i - 1 exp(-Cr)](-l) fc '(Ci?) fc '-^ 1 exp(-Ci?), 



^».n r m-r-iv- 1 1 V r(z' + » + i) y r(f" + j + i) 



c -i r(i + 1)r( |/ - i + 1) ^ T(j + l)T(l" - j + 1) 2*+'+i 



x ( k - k '- %+1 (-l) i+k ' ((R) k '~ j - 1 e?q>{-(R) 

x[(k-k' -%- l)r k - k '- l - 2 ( G xp((r){-lf +k - k '-' 1 - 1 exp(-Cr)) 
+ (Cr) fe - fe '- i - 1 (Cexp(Cr) + (-l) r+fe - fe '- l Cexp(-Cr))] , 

where T is the Gamma function 

r(n + l) = »!, 

and aJJfip) is defined by 



d V" 1 /l d N 



dx / \ x dx 



■ cxp(z) = E a k(P) xk xP ~ n ~ l CX P(^)- 
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Figure 3: Discretization of domain showing overlay of atomic center grids and 
background grid. A uniform mesh is used in the background grid while a loga- 
rithmic mesh is used inside the spheres. 



3 Numerical approximation and implementation 

In this section, we consider the discretizations for the real-space LASTO method, 
its numerical approximations, and its implementation. 

The Schrodinger equation is defined on an infinite domain and the charge 
density is rapidly decaying and smooth at large distances from the atoms. To 
accomplish an efficient discretization of the infinite domain, we consider a large 
finite adaptive domain which has fine meshes near the atoms and coarse meshes 
at large distances from them. At the edge of this finite domain, we impose 
Dirichlet boundary conditions. 

Because the basis functions are defined separately inside of the spheres (muf- 
fin tin) and outside of the spheres (interstitial region) , we have to handle differ- 
ently the regions interior and exterior to the spheres and match the solutions on 
the sphere boundaries. We use overset grids which consist of regular cubic grid 
meshes in whole domain and exponential radial grid meshes inside the spheres 
(Fig. ©. 

Inside the spheres, we represent the electronic charge density and the poten- 
tial with linear combinations of real spherical harmonics, i.e., 

p( x ) = ^2pLM(n)yLAi(h), V(x) = ^V LM {r l )yLM{f i ). (3.1) 

LM LM 

For computation, we have to restrict the A in (2.9) and the maximum for L 
in (2.14). We choose 8 as a maximum for A in (2.9) and the maximum 4 for L 
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in (3.1). These restrictions affect the accuracy of the computations which are 
also affected by the size of the finite domain, the mesh size of the regular grids 
and the radial grids. 

We next consider the implementation of the method. This consists of four 
steps: the computation of the potential, matrix generation, the solution of the 
eigenvalue problem and updating the charge density. 

To obtain the potential, we solve the Poisson equation (2.6) and (2.7) for 
the Coulomb potential in the interstitial region and inside the spheres. Because 
of the definition of Coulomb potential as the solution of the Poisson equation 
which has a source which includes delta functions at the center of the atoms, we 
use a pseudo-density p to get the Coulomb potential in the interstitial region at 
each grid point. We solve the Poisson equation 

V 2 V c (x) = 4tt(p(x)), (3.2) 
Vc{x) = 0, on the boundary (3-3) 



Ri pRi 

2 J / / „(„\ ry C _ \ Jl , 



with 

/ p{r)r^dr = q= I (p(r) - Zi8 Ci )r 2 dr 
Jo Jo 

p(Ri) = p{Ri), 

The pseudo-density p has the same zeroth multipole moments as p as boundary 
values and the same derivative also on the boundary of the spheres. 

We plot the real density and pseudo-density in radial coordinates for a single 
palladium atom in Fig. [4]and the pseudo-density and the solution of the Poisson 
equation for a palladium dimer in Fig. [5l 

To solve the Poisson equation (|3.2p . we use a finite difference scheme and 
a multigrid method [TH Q31 HH [H] , which is a well-known, rapid, and scalable 
solver of elliptic partial differential equations. 

Then, we solve the spherical Poisson equation 



d 2 2d__ L{L + 1) 
dr 2 r dr r 2 



Vt M {n) - 4^ ( P LM{n) - 8 w Zi) , (3.4) 



inside a sphere with Dirichlct boundary conditions which come from the solution 
of (HL2), i.e., 

V c {x) =Y J V£ M {r i )y LM {r i ) (3.5) 



LM 



where 



v£ M (n) = M ( £)' + [ 



Zi (I 1 



V47r V r * Ri) R 



r, 



L r Ri r Rt 



ii 
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Figure 4: Full density p and pseudodensity p in units of e/ttg versus distance in 
units of ao for Pd atom. 




(a) Pseudo-density (b) Coulomb potential 

Figure 5: The pseudo-density and Coulomb potential for a palladium dimer on 
the z = plane 



and 

V? M (Ri)= f V c {x)y LM {x)dS (3.6) 

JdSi 

pLM{Ri)= [ p{x)y LM (x)dS. (3.7) 

JdSi 

Next, we consider the generation of the overlap matrix S and the Hamilto- 
nian matrix H in the interstitial region and the spheres. From (2.1) we obtain 

(Pi(x) |-^V 2 + V(x)| ^.(x)dx = J ipi(x)(p k (x)dx. (3.8) 

Using Green's theorem the left-hand side (the Hamiltonian matrix) of (|3.8[) can 
be rewritten 

J iPiix) j-^V 2 + V(x)J p fc (x)dx 

= / ( ^,(x)V^ t (x) + (pi(x)V(x)(pk{x) ) dx - i / V^^^dx 
7n V 2 / 2 ./an « n 

Thus, the Hamiltonian matrix is symmetric with zero Dirichlet or Neumann 
boundary conditions, but not in general. Thus we consider the following sym- 
metric form derived from Green's theorem and zero Dirichlet or Neumann 
boundary condition, 

iv^i(x)V^ fe (x) + (pi(x)V(x)(p k (x)j dx = J (pi(x)ip k (x)dx. (3.9) 

The evaluation of (|3.9p requires considerable computational effort, so we use 
the following symmetrized form 

- f ^ t (x){-iv 2 U fc (x) + ^ fe (x){-iv 2 U 4 (x)dx 

2Jn 1 2 J f (3 ' 10) 

+ / ipi(x)V(x)ip k (x)dx = / (pi(x)ip k (x)dx, 



In the interstitial region, we compute the contribution of the matrices S* 
and H* by 



H, 



inlm.i'n'l'm' 



/ 4>inlm (x) < 


(-iv + v W } 







'V n' I'm 



-\v 2 + V(x) 



(3.11) 



4^inlm (x)0i / n'/'m / (x)dx, (3.12) 

In-UiSi 
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where 



</>mim(x) = r™ 1 exp(-Cr i )y; TO (f i ). 



(3.13) 



To compute the contribution of the matrices S* and H*, we use the grid or 
mid-point of each cube, i.e., 



9 1 

inlm.i n I m 



^2<t>inl m (P)<t>i'n'l'm'{Cp)Vol(C P ), 



or 



^inlm.i' n' I'm' ~ 7. / 4>inlm ( x )0i'n'i'm' (x)dx 

= (j>inlm{C)4>i'n'l'm'{C)Yo\{Cc) 



where P's are the points of a uniform grid, C"s are the center of cube in uniform 
grid, and Vol(Cp) is the volume of a cube centered P. In each case, we have to 
compute ipiimn on all grid points for each i, n, I, m. 
For the spherical coordinate operator 



d 2 2d_ 1_ 

dr 2 r dr r 2 



1 d 2 cos6 d d 2 



sin 2 6 89 2 sin < 



and 



we have 



1 d 2 coach d d 2 



sin 2 (A <96> 2 sin </> 5</> d(/> 2 



2/Jm = 0, 



--v 2 + y(x) 



[(/(J + 1) - n(n - l))r?- 8 + 2Cnr?" 2 - C'r™" 1 ] expKn)^^) 



+ V^(x)r" 1 exp(-Cr i )y im (f i ) = inim (x). 



So, we have to compute 



(3.14) 



H 1 



ImA 1 ' n'V m' 



9 1 

inlm.i' n V m' 



2 H (<?W(P)<^«'i'm'(P) + <l>i' n'l'm' (P)^nlm(P) 
P 

J2^nl m (P)(t>i'n'l' m '{P)yol(P), 



i.e., 4>inim(P) and 4>i n im(P) at each regular cubic grid point which does not 
include any spheres. 
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For the contribution of S and H inside spheres, we can compute the elements 
of the matrices with given /3's and a's as in pQ. 

S inl m A'l'm'n' = / fanlm{x)<l>i'l>m>n>(x)dX 



'Vm'jXit < 9jx\ijjX >)) 

i /■ r . ( i. 



H Ltm,i'l'm'n' = j / ( <^™( X ) 1 + ^ ^'''"'"'W 



'i'n'l'm' 



-iv 2 + l/(x) 



~t~ ~ ^ ^ {Pinlm,jX^,^i' n'l' an' ,j Xfj, &inlm,jX[i0i' ' n' 'V ' m' 



3 inlm, i' I'm' n' ' 

+ ZJ / [AnimjA M 5jA( r ) + a mimjA M 5jA( r )]^LA/( r ) 

Xn,X'n',LM,L>0 

■ \J3i'n'i>m',j\'n'9j\'(r) + a i > n n> m >,j\>i J .>gjx>(r)]r 2 drI R (\(i,\'iJ,',LM), 

where Jr(A/i, A'/i', LM) is a real Gaunt Integral. 

To solve the eigenvalue problems Htpk = CkSifk where H and S are real 
symmetric and dense matrices and S is positive definite with dimension n = 
•^basis) w e use ScaLAPACK [IB], which is well developed and optimized for 
parallel machines. 

Here, we consider updating the charge density from the solutions of Hip^ = 

We have to compute p ou t = J2e k <E F ■ IVkl" 2 + pcore m the whole region 
including both the spheres and the interstitial region. We write ipk 

inlm 

In the interstitial region, we compute 

Pout(x) = IVfc(x)| 2 , 
£ fc< B Fermi 

where 

V? fe (x) = ^ ( Pk,inim(ri~ 1 e~ Cr, yi m (f i )) 

inlm 

at all regular cubic grid points x. 
Inside the spheres, we compute 

p(x) = ^ PjLMVLMjrj) + pcore = ^ |^fe (x) | 2 + Pcorc- (3.15) 
jLM ^< B Fermi 
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From the orthogonality of g and <?, we have 



PjLM = / ip k y L M{fj)ds 

= ^ ^ ^ l Pk,inlm l PkA'n'l'vi'[PinlmJ\^9j\( r j) + ttmim. j A^5j A^i ( r j )] 

inlm i'n'l'm' jXfi j\' fi' 

[fii>n>i'm',]\>n>9]\>{rj) + Oi'n'l'm' n' 9jX' i^'(rj)]lR(LM, A^i, A'//)- 

4 Results 

To test the method we have calculated the total energy for several palladium 
clusters. 

The total energy is computed in the usual way 

OCC „ ^ 2 r 

B=-V / (fk{x.) — ip k (x.)dx+ / y Ex t(x)p(x)dx 

k JQ 1 JQ, 



2 7 n 7n |x-x'| 

where the terms are respectively the non-interacting kinetic energy, the external 
potential, the Hartree and the exchange and correlation energies. This formula 
can be simplified using (|2.ip to yield 



E = ^2e k -- / VHartrcc(x)p(x)dx- / V xc (x)p(x)dx + E xc , (4.1) 
z Jq Jn 



k 



where 



E xc = / e xc (x)p(x)dx. 



There have been many calculations of the energetics of Pd clusters in recent 
years p~8j[T9l[20l[2Tl[22]. Most calculations show that small clusters are magnetic 
(spin polarized) and for example that Pd\z has an icosohcdral or possibly a 
buckled biplanar structure [19] . Since our purpose is to test the method, we 
have compared results for non-magnetic structures in simpler geometries. 

The LASTO basis set was optimized for crystalline palladium in the face- 
centered cubic structure. It consistes of two s, two p, and two d functions plus 
one / function per atom. The £ values are given in Table [1] 

The Is through 4p states were treated as core functions. They are solved 
using the full Dirac equation in the spherical part of the potential within the 
spheres. A logarithmic radial mesh was used inside the spheres, with 429 points. 
The (uniform) mesh in the interstitial region had a spacing of 0.05 ao- Calcu- 
lations with NWchem used the SBKJC basis with an effective core potential 
given in [17] . In both sets of calculations states near the highest occupied level 
were smeared with a Gaussian of ~ O.OOlau. 
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Tabic 1: LASTO Basis and C values 



function 


C(0 


4d 


2.20 


5d 


2.60 


5s 


1.00 


5p 


1.00 


6s 


2.00 


6p 


2.00 


4f 


1.60 




Figure 6: Binding energy of Pd2 in units of Hatrees/atom as a function of the 
separation distance r in atomic units. The smooth curves are fits to a Morse 
potential. 



In Fig. Owe plot the binding energy of a palladium dimer as a function of r 
with LSD and GGA expressions for exchange and correlation energy. Here we 
computed the binding energy per atom by 

E _ E(Pd n ) - nE{Pd) 

and fitted the data with Morse potential energy function, i.e., 
E(r) = D ^ e -Mr-Ro) _ 2e ~a(r-R ) j 

where D is the minimum energy, Rq is the equilibrium internuclear distance. 

In Table [2] we compare the results with experiment [23] [24] [25] and with 
numerical results of Zhang, Ge, and Wang [20], which use plane wave basis sets 
and GGA for exchange and correlation energy, and DMol 3 [37] program which 
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Table 2: Binding energy for a Palladium (Pd) dimer with comparing previous 
results 





binding energy 


i?o 


Experiment [23 


0.37 ± 0.13 cV/atom 


4.71 au 


Experiment 1241 


0.52 ± 0.08 cV/atom 


4.71 au 


Experiment 1251 


0.698 eV/atom 


4.71 au 


LASTO (PBE-GGA) 


0.559 eV/atom 


4.802 au 


LASTO(LSD) 


0.882 eV/atom 


4.604 au 


Plane wave, GGA [20] 


0.63 eV/atom 


4.70 au 


DMoF (PBE-GGA - LCAO) [27] 


0.415 eV/atom 


4.70 au 


DMoF(PW91-GGA - LCAO) [27] 


0.538 cV/atom 


4.70 au 


DMol 3 (PW-LSD - LCAO) [27] 


0.770 cV/atom 


4.70 au 


NWChcm (PBE-GGA) [26] 


0.468 cV/atom 


4.802 au 



4.802 au 



4.72 au 



4.88 au 



Figure 7: Structures of dimer and 4- and 8- Pd clusters used to test the LASTO 
code 



uses LCAO (Linear combination of atomic orbital) with several methods for the 
exchange and correlation energy. Our results for both the internuclear distance 
and the binding energy are in agreement with previous calculations. 

In Fig. [Jj we present the structures of Pd2, Pd4, and Pds- We calculate the 
binding energy for Pd2, for Pd4, and for Pds with LSD and GGA methods for 
the exchange and correlation energy. The bond lengths were chosen to be the 
same as [20] and were not varied. In Tabled we present the results and compare 
with the results of [20J and with [25]. Again the energies show agreement with 
previous work. 



5 Conclusion 

In summary we have shown that a real space version of the linear augmented 
Slater type orbital method provides good agreement with other calculations for 
small clusters. The full treatment of core orbitals while maintaining a small 
basis set will be useful for treating heavy elements such as Pt or Au which are 
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Tabic 3: Binding energy in eV/atom for Pd2, Pd.4, and Pds 





Pd 2 


Pd 4 


Pd 8 


LASTO (LSD) 


0.861 


1.633 


2.700 


LASTO (PBE-GGA) 


0.559 


1.265 


2.189 


Plane wave (GGA) 


.473 


1.234 


1.995 


Plane wave (GGA) [2DJ 


0.63 


1.46 


1.91 


NWChcm (PBE-GGA) 


0.468 


1.405 


1.869 



important nanoscalc systems. 
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